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INTRODUCTION 

Kirchhofl' migration is an integral technique that back- 
propagates a wave field. This procedure, when applied to an 
observed scattered wave field, focuses the field at the locations 
in the medium from which it was scattered, thus producing an 
image of the medium. Among those who introduced the 
Kirchhoff integral to exploration geophysics for imaging the 
earth were Gardner et al. (1974), French (1974, 1975), and 
Schneider (1978), The technique was developed for application 
to surface seismic recording geometries and constant-velocity 
media and, subsequently, was adapted to quasi-layered veloci- 
ty structures by sequentially applying constant-velocity migra- 
tion for each layer and back-propagating the scattered field 



downward to the top of the next layer (Berry hill, 1979). This 
procedure is known as recursive migration (e.g., Berkhout, 
1982). 

In general, for two-dimensional (2-D) velocity variations, 
one cannot apply a recursive migration using the constant- 
velocity Green's functions; a nonrecursive migration must be 
employed with the Green's functions determined numerically, 
for example, by ray tracing. One recording geometry that re- 
quires the use of numerically computed Green's functions is 
the Vertical Seismic Profile (VSP). Since velocity varies with 
depth, the velocity changes across the receiver array in a 
manner equivalent to lateral velocity variation for surface seis- 
mic experiments. Such nonrecursive implementations of 
Kirchhoff migration are computationally expensive except for 
special cases. Carter and Frazer (1984) developed a fast, non- 
recursive Kirchhoff integral technique based upon Fermat's 
principle for the case of smooth lateral velocity variations. 

The goal of this paper is to describe a fast, nonrecursive 
prestack Kirchhoff migration for variable velocity media with 
large lateral velocity variations, applicable for 2-D or 2.5-D 
acoustic or elastic media, by computing the Green*s functions 
with the paraxial ray method. First, we give a brief review of 
Kirchhoff migration which describes how to back-propagate 
the wave field and form an image of the medium. Second, we 
follow with a brief description of the paraxial ray method. 
Third, we use the concept of reciprocity to implement Kirch- 
hoff migration with variable velocity using the paraxial ray 
method. The method is illustrated with synthetic VSP data 
and with VSP data collected in a borehole adjacent to a reef 
in Michigan. 



KIRCHHOFF MIGRATION 

In this section we briefiy review Kirchhoff migration and 
define an imaging condition for 2-D acoustic media. The final 
result, however, is a general form that is adaptable to 2-D or 
2.5-D acoustic or elastic media. 

The Rayleigh li integral (Berkhout, 1982) describes how 
to forward propagate the acoustic pressure field from a surface 



ABSTRACT 

A rapid nonrecursive prestack Kirchhoff migration is 
implemented (for 2-D or 2.5-D media) by computing the 
Green's functions (both traveltimes and amplitudes) in 
variable velocity media with the paraxial ray method. 
Since the paraxial ray method allows the Green's func- 
tions to be determined at points which do not lie on the 
ray, two-point ray tracing is not required. The Green's 
functions between a source or receiver location and a 
dense grid of thousands of image points can be esti- 
mated to a desired accuracy by shooting a sufficiently 
dense fan of rays. 

For a given grid of image points, the paraxial ray 
method reduces computation time by one order of mag- 
nitude compared with interpolation schemes. 

The method is illustrated using synthetic data gener- 
ated by acoustic ray tracing. Application to VSP data 
collected in a borehole adjacent to a reef in Michigan 
produces an image that clearly shows the location of the 
reef. 
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5 to a point r, 



P(r, ^) = ^ ' co)VG(r^ I r, <o) • ds. 



(1) 



where P(r, m) is the temporal Fourier transform of the pressure 
field at the point r, P(rg, co) is the temporal Fourier transform 
of the pressure field on the recording surface s at the geophone 
location , ds is a function of , G(r^ | r, to) is the Green's 
function from point r to point r^, and n is the inward unit 
vector normal to the surface. 

Migration requires propagation of the field backward into 
the medium where the field will focus at the scatterers and 
form an image of the medium. The back-propagated field can 
be expressed as 



P(r, ®) = ^ I ^frff' ©)VG*(rJr, o>) • n ds. 



(2) 



where G*(r^|r, ©) is the complex conjugate of the Green's 
function G(rg \ r, co). 

A function R(r, co) is defined as the ratio of the back- 
propagated field from the receivers P(r, co) to the incident field 
from the source P^ci^y co) (e.g., Claerbout, 1970); 



R{T, 0)) s 



P(r, (o) 



(3) 



The synthetically determined incident field Pgrcfr* ®) equals 
5(co)G(r^ I r, (o), S(co) being the source function and being the 
source location. Combining equations (2) and (3) gives the 
holographic (i.e., frequency-domain) form of Kirchhoff migra- 
tion. A pseudorefiection coefficient can be determined from 
equation (3) in the time domain by setting t = 0. The term 
"pseudo" is used because the data are not due only to nor- 
mally incident reflections. The assumption inherent in equa- 
tion (3) is that point scatterers in the medium are excited at 
the arrival time of the incident wave from the source. This is 
referred to as the ^'excitation time imaging condition" (Chang 
and McMechan, 1986). The Green's functions are expressed in 
the form 



and 



G(r, I r, CO) = ^(r, | r) exp (/(OT(r, | r)) 



G(r^ I r, ©) = Air^ \ r) exp (i03x(r^ | r)). 



(4a) 



(4b) 



where A and x are ray amplitude and traveltime, respectively. 
The subscripts s and g designate shot and geophone. For the 
far-field case, one obtains the pseudorefiection coeflicient, 



Rir) 



r)/)(r,,r,,r)ds. 



(5) 



where the weighting function W and the data D are given by 
and 

;) r 1 

(7) 



0{T,, r„ r) = ^ p|^r^. t(rjr) + T(rJr)J, 

with 9 being the angle between the normal to the surface s and 
ray (rjr). 



Note that equation (5), which was derived for 2-D acoustic 
media, has been written in a general form. The weighting func- 
tion W and the data D in equation (5) will be different, de- 
pending upon whether the imaging problem is 2-D, 2.5-D, 
acoustic, or elastic. For example, for 2-D elastic Kirchhoff 
migration, W will be a 2 x 2 matrix and D will be a two- 
component vector (Keho and Wu, 1987; Kuo and Dai, 1984). 
Also, instead of integrating along the recording surface ds, 
integration can be performed over the solid angle between the 
incoming and outgoing rays at the image point (Miller et al., 
1984). 

The weighting function W depends on the velocity struc- 
ture. The implementation of equation (5) becomes expensive 
when the Green's functions [equation (4)] in W cannot be 
determined analytically. 



THE PARAXIAL RAY METHOD 

A fast migration scheme can be obtained by computing the 
Green's functions with the paraxial ray method. The method, 
developed by Beydoun (1985a), was presented for modeling by 
Beydoun (1985b) and applied to inverse scattering problems 
by Beydoun and Keho (1986). This method is an outgrowth of 
the work on Gaussian beams by Cerveny and his colleagues 
(e.g., Cerveny et al., 1982; Cerveny et ah, 1984) and is de- 
scribed in detail by Beydoun and Keho (1987). The value of 
the paraxial ray method is its computational efficiency: the 
method does not require two-point ray tracing. Instead, the 
Green's function can be determined at any point which is 
near, in the paraxial sense, the desired ray (Beydoun and 
Keho, 1987). 

The conditions for validity of ray methods are described by 
Ben-Menahem and Beydoun (1985), and for the paraxial ray 
method by Beydoun and Keho (1987). Where dynamic ray 
tracing is not valid, for example in regions of the model where 
caustics or diffractions occur, the paraxial ray method also 
will not be valid. Therefore, the velocity model used for migra- 
tion should not contain such regions. We should also mention 
that the migration scheme we have described is not valid 
where the phase of the scattered wave field changes (e.g., due 
to postcritical reflections or caustics). This is true for nearly all 
migration methods. Therefore, we cannot benefit by using 
more sophisticated methods such as the Gaussian beam 
method for computing Green's functions in these complicated 
regions. 

The Green's function for the acoustic field [the correspond- 
ing expressions for elastic media are given by Beydoun and 
Keho (1987)] at the point M' on the ray is given by 



P{M\ (o) = Ao(M\ co) exp {ia>T(M')}. 



(8) 



The Green's function at a point M which is near but not on 
the ray can be approximated by 

P(M, CO) = qAo {M\ co) exp |i(o[^T(M') -f . (9) 

where q is an amplitude correction, which, for a point source, 
is given by 



-( 



(10) 
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5jR„. is the distance from the wavefront as determined at M' to 
the point M olT the ray, and is the Jacobian of the transfor- 
mation between Cartesian and ray-centered coordinates. The 
traveltime correction is given by the term n^KJl, where is 
the wavefront curvature at the point M' on the ray, and n is 
the perpendicular (to the ray) distance from the point M' on 
the ray to the point M near the ray. AI! of the quantities in 
equation (9) are readily available from standard dynamic ray 
tracing. 

By shooting a fan of appropriately spaced rays, Green's 
functions can be determined rapidly between the origination 
of the fan and every other point in the medium. This is ex- 
tremely valuable for migration, as well as for other inverse 
scattering problems such as variable velocity tomography or 
Born inversion [e.g., see computational requirements of Beyl- 
kin (1985) or Levy and Esmersoy (1987)] which require knowl- 
edge of Green's functions between the sources and all the 
points in the image area and between the receivers and all of 
the points in the image area. Typically, there are tens of thou- 
sands of image points and performing two-point ray tracing 
between each image point and each source or receiver would 
be prohibitively expensive. 



IMPLEMENTATION 

Implementation of Kirchhoff migration using the paraxial 
ray method was suggested by Keho (1984). Beydoun and 
Keho (1987) developed and tested the paraxial ray method for 
application to inverse scattering problems, i.e., computing 
Green's functions at a dense grid of points. Here, we imple- 
ment prestack Kirchhoff migration with the paraxial ray 
method in order to achieve an accurate and Jess compu- 
tationally expensive migration. 



Computation time for prestack Kirchhoff migration in vari- 
able velocity media can be decreased by two orders of mag- 
nitude by invoking reciprocity — replacing G(r|r^, co) with 
G(r jr, (o) and G(r|r,, co) with G(rjr, co) in equation (3), i.e., 
shooting from sources and receivers toward the image points, 
instead of shooting from image points to receivers and 
sources. Invoking reciprocity is a common procedure, e.g., see 
Hu and McMechan (1986) or Gray (1988). 

Use of the paraxial ray method further reduces computation 
time by another order of magnitude compared with com- 
puting Green's functions (traveltimes and amplitudes) from a 
2-D polynomial fitted to Green's functions determined at 
points along the rays. If only traveltimes are required, the 
reduction in computation time is less than one order of mag- 
nitude. Interpolation schemes (Hu and McMechan, 1986; 
Gray, 1988) such as polynomial fitting or spline or linear in- 
terpolation between rays will result in smooth traveltime and 
equal amplitude surfaces; however, smoothness is not an in- 
dication of accuracy. Interpolation schemes provide no infor- 
mation on how to control ray density in order to maintain 
accuracy for traveltimes and amplitudes. This lack of infor- 
mation usually results in more rays being shot than necessary. 
Unlike interpolation, the paraxial approximation is derived 
from the physics of wave propagation. Consequently, infor- 
mation is available which enables the user to control ray den- 
sity interactively (Beydoun and Keho, 1987). In most cases, a 
single ray is suflicient to determine the Green's functions at 
hundreds of image points. 

Since the data can be migrated one receiver at a time, with 
currently available field computing systems migration of VSP 
data can be performed as the data are recorded. The image 
array can be dumped at any time to look at the image as it is 
developing. Obviously, in such a situation, very little pre- 
processing can be done. In the examples that follow, however. 



Source 
0.0 1 1 — 



vi =2.5 km / s 




Offset (km) 

Fig. 2. Synthetic acoustic data generated by the paraxial ray 
Fig. 1. Acoustic ray-tracing velocity model. There are 80 method using the velocity model in Figure 1. The direct wave 
receivers on the line indicated by the five dots, and 4000 image from the source and the primary reflection from each of the 
points in the image area enclosed by the dashed rectangle. three interfaces in the image area are included. 
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very satisfactory images are obtained with minimal pre- 
processing. 

EXAMPLES 

In this section paraxial ray Kirchhoff migration (PRKM) is 
illustrated with a synthetic example and is applied to field 
data. The synthetic data example shows the results of migrat- 
ing VSP data generated from a ray-tracing model consisting of 
several dipping layers, and illustrates the effects of low ray 
density. In the field data example, a reef in Michigan is iitiaged 
using multioifset VSP data. In both of these examples the 
incident part of the wave field (first downgoing arrival) was 
removed by filtering. In the field example, ringing traces were 
killed. No other processing was applied to the data. 

Synthetic data example 

The model consists of flat and dipping interfaces (Figure 1) 
with 80 receivers and 4000 image points. The acoustic synthet- 
ic data (Figure 2) contain the incident primary reflected field 
generated using the paraxial ray method. Figure 3 shows ray 
fans from the source and from one of the receivers. The ray 
fans in this figure are very sparse and are shown only for th^ 
purpose of illustrating how each ray captures image points 
that are near it in the paraxial sense. 

Figure 4 shows the migrated data that result from using 
exact velocities and a fan of rays that are mildly sparse, i.e.. 
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Fig. 4. Migration of the synthetic data in Figure 2 using the 
exact velocities. This figure illustrates the result of migration 
with a mildly sparse fan of rays (a few percent of the image 
points are missed by each ray fan). Notice the choppy look of 
the image. Parts of the waveform in the image are zero due to 
those image points being missed by the ray fan from the 
source. This is easily corrected by shooting a denser fan of 
rays. 



Source 






Offset (km) 



Offset (km) 



P Velocity (km / s) 
0 8 



Fig. 3. Example of a fan of rays shot from the source (left) and 
from one of the 80 receivers (right). The small dots indicate 
image points. The large dots indicate the image points each 
ray captured (i.e., for which the Green's function was deter- 
mined). 




P Velocity (km t s) 
0 8 
11 I I 




Fig. 5. (a) f-wave velocity structure at the VSP site in 
Manistee County, Michigan determined from well logs and 
1-D iterative VSP waveform inversion, (b) low-pass filtered 
version of (a) used as the migration velocity. 
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OFFSET I (1220 M - OFF REEF) OFFSET H (1220 M - ON REEF) 




Fig. 6. Vertical component field data from offset I (off reef) and offset H (on reef). Notice the reflection from the reef 

area near 1400 m. 



DEPTH MIGRATION 

Offset (m) 




From Offset I From Offset H 



Fig. 7. Image created by migrating the P-to-P reflected wave field in the observed data shown in Figures 6. Notice 
strong reflector just below 1370 m on the ofTset I image. This reflector continues across to the offset H image and t 
jumps up to just above 1370 m. This indicates the location of the reef on the offset H side of the well. 
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Fig. 8. Location of source offsets I and H with respect to the proposed location of the reef as determined from 
dip-meter analysis and VSP forward modeling. The migration image (Figure 7), shown superimposed on the model, 
indicates that the reef may be somewhat closer to the borehole than indicated by the proposed model. 



such that only a few percent of the image points are missed by 
each ray fan. If an image point is missed by the ray fan from 
the source, then the source-to-image point traveltime will be 
zero. Similarly, a miss by a receiver ray fan causes the 
receiver-to-image point leg of the stacking time to be zero. 
With only a few misses, reflectors are still located accurately, 
but the image appears choppy — parts of the wavelets in the 
image are zero (e.g., see the left edge of the lowest reflector). 
This phenomenon is due to misses from the source. If the 
source-to-image point traveltime is zero, then none of the data 
from the receivers will be stacked at the correct time. Since the 
data in this example are synthetic, the traces are exactly zero 
except at the isolated times where the wavelet arrives. Thus, 
an incorrect stacking time which is too small due to a miss 
from the source will result in stacking the data at times for 
which the data are exactly zero. It is not serious if the receiver 
ray fans each miss a few percent of the image points, since 
each ray fan will miss different image points. 

As a general rule, a denser fan of rays should always be shot 
from the source than from the receivers. Since the components 



of the stacking times from the source to the image points are 
computed only once, shooting a dense fan of rays from the 
source will substantially increase the accuracy of the stacking 
time with little additional computation cost. In practice, it is 
not uncommon to perform more than one migration with the 
migration velocity updated manually before each migration. 
In such cases, we suggest that sparse fans of rays should be 
shot from the receivers until the final migration. This will 
reduce the overall computation cost, since fewer rays will be 
shot. 

Field data example 

In this example PRKM is applied to VSP data collected in 
a borehole adjacent to a reef. These data were collected in 
Manistee County, Michigan by Compagnie Generale de 
Geophysique and the Earth Resources Laboratory of the 
Massachusetts Institute of Technology as part of an experi- 
mental group shoot for reservoir delineation. The objective of 
the experiment was to determine the location of a Niagaran 
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reef with VSP data, surface seismic data, and well-log data. 
We used a subset of the VSP data for this example. The VSP 
data were recorded in a well that was about 2000 m deep with 
81 receivers separated from each other by approximately 
9,14 m and a time sampling interval of 2 ms. Two source 
locations were occupied; H was offset 1220 m from the well on 
the on-reef side of the well and I was offset 1220 m on the 
off-reef side. 

The P-wave velocity structure determined from well logs 

and l-P waveform inversion, and the macromodel used for 
migration, are shown in Figure 5. The vertical component 
data for the two source locations I and H are shown in Figure 
6. The data from each offset were migrated separately. The 
migrated images are shown juxtaposed in Figure 7. Migrating 
the data from one receiver took less than 60 s of CPU time on 
a VAX 1 1/750. In Figure 8, this image is reduced and shown 
overlain on the proposed model of the subsurface as deter- 
mined from 3-D migration of surface seismic data, seismic 
modeling, and dip-meter analysis. The VSP migration image is 
consistent with the results obtained by the other methods and 
indicates that thp reef may be somewhat closer to the borehole 
than indicated by the proposed model. 

CONCLUSIONS 

A rapid, variable velocity, prestack ICirchhoff migration is 
implemented by computing the Green's functions (traveltime 
and amplitude) with the paraxial ray method. Examples with 
synthetic and field data show that this technique produces 
high-quality images. This technique is computationally more 
efficient, by one order of magnitude, than traveltime and am- 
plitude interpolation schemes. For traveltimes only, it is faster 
by less than one order of magnitude. Since the data can be 
migrated one receiver at a time, with currently available field 
computing systems migration of VSP data can be performed 
as the data are collected. 
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FKr. 5.3-21. CMP gnthor (loft ) and the velocity sport rum (ri^lit ) as in Fignro 5.8-20. with movooni cornKiion applied to 
the data iisiiig the veloritv fujiriion thnt is posted on the spectrum. Thii? velocity function is appropriate for the nearly flat 
(events JASSoriated with thr sedimentary strata. 
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FIG. 5.3-25. (■:MF' gather (left) and the velocity spectrum (riglil) a.'s in l igmo r>.;i-23, witJi the velocity function of the 
steeply dipping events its in FiRure .'>.:j-22 jKisted on the spectninK 
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KIG. 5.3-26. CMP gather (left) ami the velocity spectriiiii (right) as In Figure 5.3-2.3. with rnoveoiil correction applied tc 
the tl;it.a using the DMO-corrected velocity fuiirtion that is posted on the spcMrtrum. 
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KIG. 5..V29. Pen urns of tlir«M' ( onunon-offset j>ecti<»ns t hat coiiiridf wit li I ho CVS p;>mHs shown in Kignre 5.3- Ui before 
DM0 roncctioji. Otfci^t.< fmn) lefl to right. 500. 1500. an<J 2500 ni. 
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FIG. r^ii-S*^. r^rtions of thrf'^ roninion-ofTsi'l soctiuns t.hnt coinridp wilti ihv CVS paiM?)$ shown in Fignrr Ti-MS aft it DKJO 
cnrn'dion. 0.!sots arc. from Kfl to riglil . 500. J 500. and 2500 m. 
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KUJ. 5.3-32. CMr^nilicr (lolt | and thn wKx it.v sjK'ft nun (right ) oi Ux niion A ns in Figure 5.3- J. S. Init after DM0 lorrection 
and ruinn)i>n-uffsrt nii«ir;uit>n. Tlir vrloriiy Innrijnjj pustt'd on the spectnini is ihf saiiio as in Figuro 5.3- 2(> derived from tht^ 
d;iia afttM* DMO rorri*rrii.Mi. 




FIG. 5.3-33. CMF giiUior (lett ) niul \ vplority sport mm (right ) at; in riRiire 5,3-H2 aflei" DMO corwtion. roinmon- offset 

iinj;r;it itui nmi iiivorsf' Ni\K) corr<:'<:iion using the vc^lority fujK tion ported on 1 lie spectrum. Tliis velocity function ii^ tin- Siunv 
ill I^igurc r>.:^-2i) ♦it?rivpjl IVouj tin- nfter DMO ( orm tioii. 
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FIC 5.C-34. CMP ^Mhvr (Irfl) and the velocity ypeclnim (right ) ,'is in Fignrf^ 5.3-33 after DMO rorroction. common-offset 
inigraticm imerst' NMQ rorn'ctiori using the velocity function postiyi on ihv spt-ctrmn in Figure .5.3-33. and finally. NMO 
' urn rtion using \hp v. lo( ity function posted on the* S))ertriini shown here and derived from tJie gather shown here after DMO 
rorreet ion and ( onnnon-offsrt inigroticm. 



Dip-Moveout Gorrc<:t.ion and Prestack Migration 




Seismic Data Analysis 




Dip-Moveout Correct iini nml Prijsiack Migration 



767 




Seismic Data Analysis 



CMP ^!;;U1ht using I ho low-vi lot iiy lunrtioij iji I'lfxiiif 
5,3-21. ovpiits I hat ('oriospo)id to tlio urntly flipping re- 
f]e( tiuns urv Hal t oned \vIut»n\.s ox e nis i])i»t con espoiKl tu 
thr s^K^eply clipj>ing reflect ionj> ,mi\ diflrni tions nre o\pr- 
corr(*cted. By apj)lying NMO correction to the CMP 
gather using the high- velocity fuviction in Figure 5.3- 
22. events that, correspond to the steepjy dipping reflec- 
tions are flattened whereas events that correspond to 
tlic gently dipping reflections and difrrfutions are tin- 
flercorrected. 

The velocity spectrum ni the saine HUcdysis locn- 
tit>n lus for Figure 5.3-20 after tlie application of DM0 
correction shows a single, tnianibignons velocity trend 
sliown in Figure 5.3-23. We have resolved tlie dip ef- 
fect on stacking velocities anil thus eliminated the mul- 
tiple number of velocity picks as in Figtire 5.3-20. .As a 
result of the partial migration effect of Dh\0 correction. 
entTgy is moved from one CMP location io another. 
Consefjuently. reflection-point .smearing is lemoveti and 
evfiits on CMP gathers imj>ro\ ise a horizontalh- laytTCii 
" model. 

l^'or comparison, the low-velocity function as.soci- 
ated with the gently dipping reflections in Figure 5.3- 
21 is superimposed on the velocity spectrum after DMO 
correction iti Figure 5.3-2-^. As anticipated, the picks im- 
plied by the velocity spectrum after DMO ( orrection are 
fairly close to the velocity function associated witl) the 
gently dipi>ing reflections. SimiJaily. the high-velocity 
function associatixl with the steeply dipping reflections 
and diffractions in Figun* 5.3-22 is siiperimposed on the 
velocity spectrum after DMO correction in Figure 5.3- 
25 In this case, the picks imt>]ied by the velocity spec- 
trmn after DMO correction are significa})tly lower titan 
the velocity picks associatt^l with tlie steeply (lipping 
r(*lIections. 

The velocity analysi.s is rej)e;ited after DMO eor- 
rer lion as shown in Figure 5.3-26. The picked velocities 
are then used to stack the data (Figure 5.3-27). Com- 
pare with the couAentiona) CMP stack in Figure 5.3-1-4 
and note that the DMO stack shown in Figtire 5.3-27 
has preserved the events with conl^ictijig <]ips with dif- 
ferent stacking velocities the gently dipping reflec- 
tions with tlie velocity functiot) posted in Figure 5.3-21 
ajid the steeply clipping fault -plane reflections with the 
velocity function posted in Fignre 5.3-22. The result- 
ing migrated section in Figure 5.3-28 shows a superior 
image of tlie fault blocks as compared to the migrated 
section without DMO correction (Fignre 5.3-15). 

For prestack time migration, we shall follow the 
same sequence as for the salt -flank <lata set. 

(a) vSort Die moveout-eon ected CMP gathers to 
connTion-oiiset sections. Figure 5.3-29 shows j)i)r- 
tions of three selected con nnon- offset sirtioiis. Note 



tlie presence of steeply dipping fViuh -plane reflec- 
tions. Without DMO correction. CMP slacking 
fails to preserve t hese reflections as shown in Figure 
5.3- J 1. 

(bl Apply DMO correction to each coin njon-olf set sec- 
tion. Fignre 5.3-30 shows portions of the three 
selected common-offset sections as in Figure 5.3- 
29 after DMO correction. Note the steeply dip- 
ping fault-plane reflections. With DMO correction. 
CMP stacking preserves these reflections as shown 
in Figure 5.3-27. 

(c) Following DMO correction, each conmion-offset 
section is assumed to b(? a rei^lica of a zero-offset 
section, and thus, can be migrat<.^d using a zero 
offset njigration algorithm. 1)) this case, a single, 
vertically varying velocity function was used lo 
migrate the NMO- and DMO- corrected common- 
ofTset sectiojis. Figure 5.3-31 shows j)or lions of the 
three selected common- offset sections as in Figure 
5.3-30 after conmion-ofTsei migration. 

(d) Sort the migrated common-offset sections into 
CMP gathers. Events on these gathers are now as- 
simiingly close to correct subsurface positions after 
migration. Figure 5.3-32 shows a CMP gather after 
common-offset migration . 

(e) Apply inverse NMO correction using the veloci- 
ties picked before DMO correction and perform ve- 
locity analysis. The velocity spectrum computed 
from migiated (lata is showti in Figure 5.3-32 with 
the DMO velocities posted on it fo) comparison. 
The CMP gather after inverse NMO correction us- 
ing the >elocity function shown in Fignre 5.3-32 is 
shown in F'igrne 5.3-33. 

(f ] Apply NMO correction using tlie velocity picks af- 
ter migration (Fignre 5.3-34). The velocity field 
deriwil from the post-migration velocity picks ;vs 
in Figure 5.3-34 is shown in Figure 5.3-35. and 
the stack basc»d on this velocity field is shown in 
Figme 5.3-36. This stack indeed is equivalent to 
prestack time-migrated section. Note that, how- 
ovcT. common-offset migration actually was done 
using a single, vertically varying \elocity function 
as jji step (c). 

fg) To obtain the migratc=»d section with tlie velocity 
field (Figure 5.3-35) derived after common-offset 
migration, first perform inverse migration of the 
resulting stack from step (f) using the SJinie ve- 
locity fnnclion as in stej) (c) to obtain a zerooffset 
section ectnivalent to an umnigrated stac^k as shown 
in Figme 5.3-37. Then, migrate this zero-offset sec- 
tion using the p».»st-migrati(m velocity field shown 
ill Figure 5.3-35 to obtain tht^ final restdt frotn 
prestack lime juigration secpience described here 
(Figm-e 5.3-38). 
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Common-Reflection-Point versus 
Commoi -TlefJection-Surface Stacking 

C.'ojivcMHiunal stacking is basnfl on the notion of i\ roni- 
moii midpoint (CMP) and migration is based on tho no- 
tion of 9 connnon reflection point (CRP). In both en.ses. 
we rtssuijie that reflections are represented bv hyperbo 
ln,s and reflectors are represented by points. Consider a 
seismic line from tlje Alberta Plains of Western Canada. 
The snbsnrfare is jnst as flat as the surface over which 
von iinv(- record<»d the data. When stacking the data. 
>ou can ahnost picture smnmiiig of the ajnp]itudt?s in 
a CMl^ "ather over all offsets along a Jiyperbolic tra- 
jectory cussoriated with a zero-offset time and the re- 
sulting stneked anjplitude being placed at a poivt n- 
flcclor whm the CMP raypaths converge conveniently. 
Wlien migrating the stacked data, again, you sum the 
aniplitmles along a hyperbolic trajectory and place the 
result ;n the apex of the hyperbola. You conveniently 
associrtfe the apex of the latter hyperbola with a poini 
: 'r.-. or situated on the reflector. Whether it is a n> 
iiu iiun hyperbola associated with a point reflector or a 
(ii^rrariioii In perbola associated with a point diffractor. 
the proc(>ss of stacking and migrating the data involves 
sunnn.-tioii of amplitudes along hyperbolas and placing 
the resulting stun to a point in tiie subsurface. 

Now consider ;i seismic line from the Alberta 
Foothills where the Ciuscaded flanks of the Rocky Motm- 
tains rise steej^Iy. TJie subsurface is just as steep as 
the surface over which you have recorded the data. A 
diffractor .still is a point whether it is in the Plains 
»r -n Ml * Foothills, and you can still think of diffrac- 
.iojis as hyperbolas so loJig ni? ihey are situated bek)w 
a simple t»\erburden. Fortunately, yon can also think of 
reflections as hyperbolas. No longer. howevei\ can yoti 
associate a reflectiou hyperbola on your CjMP gather 
with a single pt)int reflector: instead, yon have reflec- 
tion points dispersed along the reflector (Sect ion E }') 
This is when you ha\e to introdtice DM0 correction to 
account for the reflection-point dispersal (Section 5.1). 
Once the reflection-point dispersal is removed, the re- 
sulting stack can be considered equivalent to a zero- 
oifsei .section winch you can migrate, again, tising the 
liyperl>olic sunnnation rule. Thus yon have been able 
t(. overcome the Foothills problem of steeply dipping 
events. 

Finally, consider a seisinic line from the Canadiaii 
I'hrust Belt west, of (he Foothills. The subsurlace is 
just ns I onipl<\N ;is it appears on t he smface over which 
you liMve recordefl th»^ «l»trt. A diflVnctor situated below 
the conipliW overbiu^len srrnct\u»> cau.sed by the over- 

Ihnist itritHiirs .stays as a poini: no longor. however. 



ran yon associate it with a hyj)erboljc traveltime. In- 
stead, when migrating the data, you have to deal with 
a lomplex. distort tnl traveltime trajectory . And neither 
can you associate the refle(*tiojis with hyperbolic trav- 
eltimes. Instead. wh(Mi stackijig the data, you have to 
deal witli a n on hyperbolic moveout trajectory associ- 
ated with many reflection points scattered around in 
the subsurface (de Bazelaire. 1988). 

So the simple hvperbolic imd poijjt rules of the 
PI ains or the Foothills are no longer applicable in the 
Thrust Belt. Toovercojne the first problem • migration 
of data in the presence of strong lateral velocity varia- 
tions a.ssociated with complex overburdc.'n structures in 
the TInust Belt. >ou may decide to do the imaging in 
depth (Section 8.0) i)istead of imaging in time (Section 
1.0). Earth imaging in depth reqnires eaiih modeling 
in dej)th (Section 0.0) — a challenge nnich higher than 
the ijnaging itself. To overcome flie second problem • 
stacking of data with nonhyperbolic moveout. you may 
eombine it with ihv first problen) and pursue a rigorous 
solution l)y doing the imaging not jnst in depth bnt also 
before stack. 

Note that, insofar as stacking and imaging in time 
or in dc^pth. we choose to map events to common- 
reflection points. By way of DM0 correction, we map 
tnents to common- reflect ion points in their unmigraied 
posifiorts. Similarly, by way of prestack time migration, 
we m;tp events to connnon-reflection points in their rni- 
urnial positions. De Bazelaire (1988) and Gelchinsky 
(1988) pointed out tliat. rather than associating the 
recorded data with coimnon'refiecliov poirih. we sliould 
associate tlie reconJeci data wit h (:onimo7)-rcfleciion sur- 
fncefi. As such, when stacking the data, rather thaji con- 
stricting the summatioj) of amplitudes to within a sin- 
gle CMP gather, data may he focused to a common- 
reflection surface (CHS) using nndtiple shots and re- 
ceivers (Gelchinsky. 1988: and IJocht. 1998). 

The theory and practice of the CRS stacking 
method (Gelchinsky. 198cS: Tygel et al.. 1997: Gelchin- 
sky c\ al.. 1999a.b: Hocht. 1998: Berkovitch et al.. J998: 
Landa et a!.. 1999) suggest that it is indec^d related to 
piestack time migration when the medium velocities are 
judged to be within the acceptable bounds for the lat- 
ter. The CBP gatliers from i)restack time migration and 
the CRS gathers both are created by including in the 
summation aperture multiple lunnbers of CMP gath- 
ers. The differeuie is that the CRP gathers are associ- 
ated wit h events in their rnigraled pos^iiions, whereas the 
CRS gathers «re associate* I with events in their unmi- 
gralcd pofiiiioit.s. So. niigmtior of tlie CRS stark .shoiiJd 

resemble I he section from prestack time migratio]i. 
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FIG. 5.3-39. A rDninion-roliort ion-surf j»ro (CliSi stark. (DatJi roiirtr^y Sfiii<li Aminro.) 
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FIG. 5.3-40. Minrafion Dl tlie roinmon-roflociioii-siirlarp (CRS) stark shown in Figure r)..^39. 
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FIC. 5.3-41. A roinnK,n-i.fhvlion-,>o»nt (CRP) .tark of .he ,lai:. i,, n^.m- o.:^.|0 ,|.,ivr.i f><M„ pivslnvk tiin. n.igralion. 



DjI>-Movec>ut Correction nnrl Prestack Migration 



773 




CO 



3 

to 
IZ 

.E 

c 

is 
c 

"to 

X 

CO 



CO 



c 
S 

c 

c 



U4 



771 



Seismic Data Analysis 



0.00 



o.so 



10001 
_l 



c rp_gather 
0.00 



1.00-- 




FIG. 5.3-43. A r.«n.non-rotl.H-..ion.p<,i..t (('RP) gather asso. iate,! with th. CHP stack .shown in Figure 5.3-41. 
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FIG. r».4-l. (a) A CMP gather containing a single rotler- 
tur ill a conslani -velocity niedinni: (b) velocity spectnun 
(lerivttl by niigratitiR tlie CMP gather using a number of 
constniit velorilies and display ing ihe wro-ofTset trace from 
' •M-h migration; (r) velocity spectnmi dorivtHl by NMO cor- 
roding ajid stacking tlie CMP gather usin^ the same range 
t,f constant velocities us in (b). Eac h trace in (b) is a zero- 
offset trace, while each troce iji (c) is a CMP stacked trace. 

Figure 5,;V39 shows a CRS slack from an area with 
low-relief St ntctures. The field data were recorded with 
2 !0-fold coverage. I^or CRS processing, the rjegative off- 
sei.-^ vere discarded and the near one-half of tlie positive 
offset.s were kept: thtus the redticed data used to create 
Ihe Ci^S stack had 6()-foUl coverage. For prestack time 
luigraviun. the original 2 JO-fold data were first reduced 
to ()()-fold by partial stacking. The jnigrated CRS stack 
shov n in Figure f).3-10 is comparable to the CRP slack 
(\ . rd frfnij presiack time migration show^i in Figure 
5.3- J L In the present example, five CMP gatliers were 
used to create one CRS gather willi 300 traces as shown 
in iMgun? r).;V'42. Tlie 60-fo]<l CRP gather at the same 
location as for the CRS gat Iter of Figure 5.3-12 is shown 
in Figure 5.3- J3. 

Practical rehnejuents to various approaches to CRS 
stackijjg and imaging, including extensiojis to 3-D. are 
mider current investigation. As to what extent the tetli- 
niqne will l)e applied loutinely to seismic data remains 
to be i'i rn. 



5.4 MIGRATION VELOCITY ANALYSIS 

Velocity estimation. CMP stacking, and migration gen- 
erally are Cijnsidered independent pvoce.sses. However, 
they all have a comnion theoretical Ix^e — the scalar 
wave equation. Solution of this equation allows dowjj- 
ward extrapolation of a seismic wavefield recorded at 
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FJG.5.'l-2. Migrations of the CMP gather in Figure h.i- 
la using velocities that are (a) lower than the medium, (b) 
the same as the mpdium. and (c) higher than the medium 
velocity. (Note the wraparound effect in (c) that la inherent 
to the phase-shift method usoti in this analysis.) 

the eartli's smTace. in turn, downward extrapolation 
provides a ba^is for CMP stacking and migration (Clay- 
ton. 1978: Vihna/ and Claerbout. 1980). Because the 
processes of CMP stacking and migration require ve- 
locity infornjatitni. they also can be used to obtain a 
velocity estimate (Taner and Koehler. 1969: Gardner et 
al.. 197 I). 

For example, consider the problem of conventional 
vek>city estimation for stacking. Figure 5.4-1 a sliows a 
CMP gatlier over a single horizontal reflector. Select a 
constant velocity, apply NMO correction, and stack the 
traces in the gather. Next, place this stack trace onto 
the plane of velo{nty versus two-way zero-offset vertical 
time (r.r) ;ls shown in Figme 5.1-lc. By stacking tJie 
gather with different constant velocity values, the entire 
(r.r) plane is filled with stacked amplitudes. (In this 
section. th<* variable r is used specifically to define the 
vertical twoway zero-olfset time, r ^ 2z/i\) 

We now consider the migration process. For a hor- 
izontally stratified earth, as in Figure 5.4-la. we cannot 
distinguish betw'cen a CMP gather and a common-shot 
gather (CSG). Moreover, since a CSG is a true wavefield 
createtl by a single shot and recorded by many receivers, 
it seems reasonable that the CMP gather in Figure 
5.4- la can be jnigrated by treating the reflection hy- 
perbola as if it were a diffraction hyperbola. Assuming 
there is no velocity information, we migrate with various 
trial constajit velocities and evaluate the restilts. Figure 
.^>.4-2 shows three different attenjpts at migrating the 
CMP gather of Figure 5.4-la. In one attempt (Figure 
5.4-2a). too low a velocity w{L<^ used; hence., the event 
was undertnigrated. In another attempt, too high a ve- 
locity was used and the event was overmigrated (Figure 
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so^>rMARy 

A sigiuficAiit dcgra^lation iii the <iuality of Kirchhoff mi- 
gration images can arise when I he migralion operator sum- 
iTitiUdn trajectory is too «tccp for a given input ftRismic trace 
spading and frequency content. This phcnotuciiou is callai 
"^operator abasing" and Is strictly distinct from and indepen- 
dent of data aliasing, Unfortuiiatciy, tlio conditions winch 
lead to significant operator aliasing are common iu raosl 3 D 
KtrchhofTimgratiou work, due to tlic sparfif* and irregular Jia- 
ture of 3'D acquisition goomctrifjs. am] the need to mtnimize 
the siae of 3-D migration input data volumes. 

Wc aali-ajias the niigration op<'.ratars "oa tho fly** 
N-|X)int triangle iilter.^. Usinf; a Z-transform representation, 

rediifie each K point triangle filter operation to an efttcient 
thr«45-potnt filter operation, with the added nunor overhead 
of a causal and auticausal Intogratlou of the input trace dato 
prior to migration. Our anti-alia^ed migration is compnta- 
tioiiaily eHlcicnl since it requires about ifie .same number of 
lloating point operations as a conventional Kirchhoff inigra- 
ilon» and is eonscmtiv« in its mcmorj' use since it docA not 
require storage of multiple trace- daUt copies to impicmcjtt 
the anti-ah'asing. Our anti-aliased migration is accurate in 
that it applies local aiiti-aliafi filters sei)arate1y for each of 
tiic many migiauon operators which pass through a^giycu 
input data point. 

Wc demonstrate our method with a 3'D implemcutaliou 
on a maasivdy parallel Connection Maciiinc CM'5, and com* 
paxe our new ontl^aliascd KirchhofT migration to a convcu* 
tional. aperturc*wcighted Kirclthoir migration In an applica- 
tion to a 3^D marine seismic data example. Our results indi* 
cate that our anti-aliasing method greatly enhances the 3-D 
resoUiiion of steep salt*sedimcnt intei faces and faults, and 
suppres^jes false reflections caused by convcntioa<il Kiichhoflf- 
migration aliasing artifacts. Our method is applicable to 
Other space*ttme I^irchhoif operations such as NMO velocity- 
stacking and DM0. 

INTRODUCTION 

The process of seismic migration invdvi.'s m niap[Mug iif 
the input seismic trace data by a mi/yration operator lo (he 
outjjul stfuctund iu»age. SjKUial .ili^w>!ng of the migralion 
process C4UI be probieinaiic in all three domains: data, op- 
erator aiid image. These three types of spatial aliasing are 
distinct and independent, although i\\ty are often confused. 



Aliasing dcfiuGd 

DatA aliasing occurs when there arc seismic frequencies / 
and dips 0 along tcAcction and dilTraction events iu recorded 
seismic data whidi violate the following criterion: J < Jd^ 
Vrj{Af»\nOrd^x). The term U is the maximum unaliascd fre- 
cjuency in the data for a given trace spacing Az» a veloc- 
ity Vr at the receiver location, and a local plane- wave inci- 
dent angle (?y at the recording surface. Once data have been 
aliased during the recording process, it isdiflicult to suppress 
subsequent processing artifacts except to rcshoot the survey 
at a finer group spacing, or try to interpolate the original 
data without interpolating the aliases, e.g., Spit2 (1991) and 
Claerbotit (1992). Vilmaz (1988) ipvcs a good over\*iew of 
data diasiiig and its effects on all methods of inigratiun. 

imoffc aliaatng occurs when the output sampling of the 
image space is too coarse to properly represent migrated dips, 
I'his is fiKistly a coNiuotir prohlc^in and can be rumedicd by 
choosing aii appropriate image mesh upon which to estimate 
the migrated structural image. Note that the Kirchholf- 
migrated values at each mesh poitit in an aliased image spaco 
arc uncontaminatcd and correct (assuming no daia or opera- 
tor aliasing). For example, a single migratefi output twee at 
a well location may be entirely accurate without having toi.-S' 
timalfi any adjacent migrated traces, litis "target-oriented" 
capabUity of Kirch hoff migration is not p<»sible with spec- 
tral or finite-difTerence migration methods. Kirddioff image 
aliasing is only a problem from an aesthetic point of vkw, 
and wJicn the migrated inaa^ is needed as input to suIjsc- 
quent multi-diannel processing. 

Opcrat<ir ciiasing occurs when the operator dip along the 
migration summation tfajector>- is too steep ifor a^ven ini»ut 
seismic trace spacing and frequency conteiit. Spatial aliasing 
of KirchhofT iiitcgrol operators is likely because the operators 
arc designed in the discrete space-thnc domain rather than 
tlic sj>cclral <loinain. Groat care must be taken lo ensure that 
the f-k spectrum of a giv^n space-time operator contains no 
aliased cnerg>' wrapf»cd l>eyo«d IJie Nyquist frequency and 
wavenumbers. This implies that space*time operators which 
are spatially nndersampletl and contain significant st««'p-dtp 
and high-fre<tuency energy an.- likely to suircf from severe 
operator aliasing. In contrast, migration operators designed 
explicitly in the f-k domain are unaliased by deHnition, and 
operators designed in the /•;f domain arc easily modified to 
suppress spatial aliasing. 

Previous work 

Convcntiotially, anti-aliasing KirchhoJf operators can 
bi! partially accomplislH^I by tian: int«?rjKjlation, <-»per(ure 
control, ajid operator dip fdtering. Trace inlcrpoJalion dc- 
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creases the spatial ^ampie ijiicrval .\nd thus increases the 
spatial Nyquist frcijuency, but at the pottiitially pxoltibiiivc 
cost of increasing the total 3-D nugration inpwt load {e.g.. 
Yilmaz, JfiflS). Migration aperture control and operator dip 
fiUcring have the effect of auppressiiig aliased steep-dip con- 
tribiilioiis to the output imrige. UnfortMiiatdy, xUv.sc s«)>- 
jiresscd steep dips may be precisely the eveiits of iiUertrjit 
wlien imagifig near-vertical structure such a^ salt flanks and 
faults. 

Haie (1991) proposed a iaethod to anti-alias Ktrclihotf 
dip-inoveout (DiVIO) operators by replicating time-shifted 
Versions of the input saisntic traces to approximate a suin- 
j&atioa along an cqual-tiuic-sainplud version of the operator. 
Hale's subsequent nearest- neighbor spatial interpolation im- 
plies a boxcar averaging of the input trace values along the 
lirne direction, wliereas we use triangk filtcriag to increase 
accuracy at steep migration dip^. Hale's test of summing 
DMO impulse response:* and matching the result to ihe ong- 
inal impulse aisures both operator aJiasin^ and image alias- 
ing are suppressed, which may invoke mor^ lowpass-Bltering 
and loss of resolution than is strictly necessary for the case 
of migration operator anti-aliasing alone. Hate's replication 
of time-shifted input traces and spatial interpolations are 
InefHcient from l>oih a floating-point a!ul memory perspec- 
tive for 3-D migration algorithms, such as ours, that process 
thousands of traces simultaneously. 

Gray (1992) proposed an approximate but cost-effective 
method to anti-alias Kirchhoff mifjratioti operators by pei- 
forming aporture-depeudent Itnvpass temporal fdlering of ihi; 
input trace data. He suggeslcd that approximately three 
copied of the input data be lowpass-fiitere^l at diflerent hijf;li- 
cut frequencies. During migration the wide-aperture por- 
tions of the migration operator would be drawn increasini^iy 
from the lowest* frequency-content filtered trace copies. This 
corresponds to the noticr. that the wide-aperture portions of 
the migration operator have the steepest dt}ys and tiierefoie 
only the lowest frequencies will be unaliased, W<» more accu- 
rately calculate the coniplete vector of high-cut frequencies 
separately for each diffraction trajectory, ajid desi^^n and ap- 
ply our anti-alias lilters on the fly accordingly. Our method 
is as cost-effective as Cray*s approcKii, but docs not require 
the memor>'-inicnsivc storage of multiple data copies to im- 
plement the anti-aliasing. 

We present an improved method for anti-aliasing Kirch- 
hoff migrations. We first rierive the relevant operator alias- 
ing criteria anil then present our ellicient method of auti- 
aliasing the migration operator by local triangle lilteriug of 
the seismic trace. Finally, we compare applications of our 
anti-aliased Kirchhoff migration to conventional aperture- 
w<i]ghuxi Kircldioff .migration on a 3-D marine seismic data 
exautple. 

AN'ri-AIJASING THEORY 

The operator anti-aliasing criterion defines the inaximum 
unaliaseij temporal frequency in the i>eismic data at the loca- 



tion that the operator trajectory intercepts a seismic trace, 
and 15 a function of the local operator dip and the trace spac- 
ing. In order to remain unaliased, the sequence of seismic 
trace samples summed along the operator trajectory must 
satisfy the dip Nyquist saiupling criterion: 



U) 



The factor AT* is the moveout time along the migration op- 
erator between a<ljacent traces. The term dti,/dp is the local 
spatial derivative of the migration operator at the point of 
intersection with the seismic trace, and A/> is the seismic 
trace spacing alon^j the recording surface. The anti-aliasing 
criterion (1) suggests local lowpass filtering of the seismic 
trac* to ensure that no frequencies larger than /^^, are in- 
correctly migrated into tlie output image. As discussed later, 
we implement the anti-aliasing iowpass fdtering cflicienlly 
and robustly with local triangle fd ters. 

KiTcIihoff time migration 

It. is instructive to evaluate the anti-aliasing cciveriou for 
tlie special case of Kirchhoff Ume mifjration, which can be 
cxtende<l to the case of Kirchhofl* dept)i migration. For time 
niigrailoit, wc derive an expression for the anti-aliasing cri- 
terion (1) by anaiyiicaliy evaluating the operator derivative 
term Otk/dp. Wc begin by considering the analytic expres- 
sion of the Kirchhoif 3-D prestack time migration operator, 
which has the form 



The tenn h is «he total two-way Kirchhoff migraiio4i re- 
netlion travchime. 7 is the one-way vertical traveSlime or 
l^seudodepth to the refieciion point, and v is the rms jiiigra- 
tion velocity that varies as a function of space and p.<$eudo- 
depth (a;,;/, r). The lengtlis and p^. are the distances mea- 
;;ured along the planar recording surface between the source 
or receiver, respectively, to the vertical projection of the im- 
age point on the recording surface, exj)rc$scd as 



tkiuation (3) assumes. that the recording surface h a hori- 
zontal plane at a constant datum depth xo. The subscript 
ir is a wildcard for either the source or receiver coordinate 
subscripts s and r, and the subscript t refers to the surface 
location of the subsurface image point pix>jccted up onto the 
recording surface. 

The spatial derivative of (2) is needed to determine the 
migration operator aliasing criterion, which can be derived 
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evaluate the cfTcctivc rms spatial sampling interval A/j 
for the prcstack migrAlioft operator as 



where viv, define <far* and c/j^, in the prestack case as 



and 



(5) 



(6) 



(7) 



The term^ Ax and A3? are the true iiiliuc and crwyjline sam- 
pling ifttervals of the sdsmic tracts data, aiid arc assumed to 
he regular liut po&siUy uawjual iji length, A more compli- 
CAted analysis is required for irregular trace spadngs, which 
this j>*iper docs not address. 

Given the Analytic expressions for the spatial derivative 
of the KirchhofT tiiiic migration operator f74/ri/), and tht: 
definition of Lbe areal trace spacing interval A/., vtfi derive 
an analytic expression for the time miyratioo opcnitor anli. 
aliasing criterion. For the prestack case, the anti-aliasinf 
criterion is 
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(8) 



where the prcstack dfrfnutkni of Afi given by (ft) « (?} must 

be liAed. 

Kirchhoff depth migration 

In contrast to the time migration method afjovc, the 
depth aiigration anti-aliasing criterion dous not have a gen- 
eral analytic solution. KirchJioft depth nugration operators 
m v(x,y^ z) luc-dia do not liave simple l)yi>erholic forms, and 
must be evaluated numerically llnough heterogeneous migra. 
tjon velocity models, Iio*%-vcr, propose three methods to 
provide satisfactory depth migration anti-aliasuig results. 

Tiio finil and mo,<t straiglitfonvard method is simply to 
use the rms velocity field v and the lime migration anti- 
aliasing criterion (8), of the previous section, to design tlic 
ami-idiasittg lowpass filters along the non-hyperbolic depth 
migration operators. Tliis approach is robust because the 
djp .-Uong tlic time inigration operator is a reasoiiablc ap- 
proximation to the averse dip Oti/Op along the tme depth 
migration trajectory. 

The s^ond depth nngratiga anij^aliasing technique ap- 
plies lo Kirchhoir depth migration driven by traveltimc ta- 
bles. In this cas<s the operator derivative Otl/dp in the gen- 
eral anti-aliasing criterion (1) can be evaluated bv differ- 
encing the source and receiver trav^ltimei^ at the trace loca- 
tion with respect to an adjacent fsonra^ ;nid receiver location. 
This approach is more accurate than the lime migration ajv 
proximation, hut requires t..xlra memory a«d= r/0 overhead 



to load the adjacent Iraveltime tables, or extra compuia- 
liofi time if the adjacent travchimcs arc intcrpolatc<i from a 
sparse set of Iraveltime tables already loa^led into memory. 

The third depth migrAtton anti-aliasing technique applies 
to Kirchhoff depth migration in which the travellimc values 
are computed by raylracing on the lly. In this case, an extra 
set of rays can be tr,xx-d from adjacent source and receiver 
locations down to all subsurface points, or interpolated from 
a single set of dynamic rays by imaxial interpolation. Then 
the difference in ray arrival times can he useri to approxi- 
mate the operator derivative dh/dp in (i). This method is 
more computationally intensi\*e thaji the travellime table ap- 
proadi, but has significantly snjaller memory rcxiuifemcnls. 

In the sections that follow, ve assunm that the opera- 
tor derivative of equation (I) is a known quantity, 
whether for time. 01 depth migration, since it can be cilcu- 
lated by any of the three methods described above. 

ANTI-ALIAS FILTER DESIGN 

Given that we know the maximum unaUaj;ed fr..Hjucney 
fmc* at each point along the migration operator, we ncejj 
to design a method to locally lowpass filter the s^^tKuiic data 
in order to satisfy the aMti.ali.v>ing cfiieriou (1). We iiuple- 
mcnl the local lowpass filtLM^ n.s triangular smootJdni', The 
'^transform representation of on arbilrai>- Appoint triangle 
filter can be expressed as ^ 



^ 2 - » « 



(9) 



where the triangle lengt;h /V 2k \- i, and the filter scaling 
cocmacal tv {k -r 1)'^. The amplitude spectrum of the 
triangle filter can b«j d<;rivcd as 



(10) 



whicit is the discrete version of a coiitinuutjs siiic^ function, 
and has spectral notcht^s at the frequencies 



(11) 



We design the triangle filler leni(1h /V such that the maxi^ 
mum unahased frequency is equal to the first notch in 
the triajigle spectrum at frequency Equating the ami* 
ahasm^ criterion (I) lo the fin?t triangle filter notch fre- 
quency of {li), wc derive the appropriate anti-alSasinK trian- 
gle filter length N to be 




The valnt.' of the. local migration 0}>eraior dip mk/r}p Is given 
by erjuation (4) in the case of time migration, or by the mclh^ 
ods discussed in the previous s^xlioii oh tlepth migration. 
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Furthermore, note that the denominator of^Cc) repre- 
sents a causal integration of the trace with 1/(1 - z), followed 
by an aiiticausal trace iiitegratioii with 1/(1 z'^^). The nu- 
merator of g{z) is a gapped three-point Laplacian operator, 
where the kngU: of the gap {k -I- 1) determines the triaa* 
gle length A' and high^cut frequency Therefore, if the 
input trace data are integrated twice before migration, an 
arbitrary' A'-point ani]^a!i<is triangle filter can be designed 
and applied on the fty for the cost of a three-point opera- 
tor. This observation olfers a large saving* in computational 
<ijTort when applying the triangle filters. 

DATA EXAMPLES 

We compare our antj-aliased KircJiholT migration to a 
conventional aperture^weightcd.Kirchhofrxnigration, both tm- 
pleiuented on a Connection Machine CM-S. Preliniinafy 3-D 
results were shown by Lumiey and Claerboiit {199.3).' Fig- 
ure 1 shor^^s how onr anti-aliasing method enhances a 3-D 
KircJihoJT migration image of a salt intrusion, given a sparse 
60m-blnned input datasct. Wc vili di&cuss 2-D and 3-D data 
cxarapics in more detail at the oral presentation, 

CONCLUSIONS 

We jiresent an anU-aUasiuttg method for Kirchhoff migra- 
tion operators, based on local tiiajigle filtering. The Ap- 
point trianglfs are emciently applied as 3-point ftlter.^ after 
causal and anticausal integration of the seismic trace data. 
This method is well suited for parallel processing algorithms, 
where storage of multijjje copies of various high-cut filtered 
versions of the input data is not often practical. We com- 
pare the anti-aliasing mctliod to a standard Kirohhoff 3-0 
migration, using a salt intrusion data set from the Gulf of 
Mexico. 



Our results indicate that careful operator anti-aliasijig can 
enhance sieep-dip rcrtections from faults and salt-s^diraent 
interfaces during the process of 3-D seismic migration imag- 
ing. Our method is applicable to other Kirchhoff space-ti me 
operators such as NMO velocity stacking and DM0. 
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Figure 1: Kirdiholf :m; migralvd tinw sUctii at 1,V44 seomdh. The kft imnel is 
a&iaiidard ICi rch holf a-.D migration with a (iO-degree migration aperture and cos'-*^ 
operator dip hitcr, the right panel is our new ajUi-aliascd migration mctbod. 
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FREQUENCY-SELECTIVE DESIGN OF THE 
KIRCHHOFF MIGRATION OPERATOR^ 

SAMUEL H. GRAY* 



Abstract 

Gkay, SJEl 1992. Frequencjy-setoctivo dAsign of the Kirctooffmignition opewtor. Geojkytleat 
Prospecting 46, 565-511. 

Integral nugration teoliniqu« perform a «Jm over a» aperture of input traces to obtain 
output at a single point The length of the aperture is limited by a SP^.'^^ 
^ch typically prohibits imagiag very steep dips at very Kghfrequenoes ^fo"* 8J«^« 
seme migratiOT artifects (migration operator aliasing). For nme^omam KfT'^lT^" 
tion, this can be a fatal shortcoming. The standard way to axJdress this probtao is to mte- 
polate traces spatiaUy before migration. This reduces the trace spacing, thereby hicreasmg Hie 
lEmuency content idndi can be migrated without aliasing at steep dips. 

An alternative remedy to the operator aliasing probtein is to modify the ptase response 
of the Kirdbhoff migralioo operator. This operator is frequency-sdective across the mxgra&OD 
apertuie: it passes all temporal ftequsndcs of the input traces in the innermost portion of Uw 
apertme (jrferring to the shallow dips), and graduaUy cute out the l»if»J»4««a« « ^ 
approaches the outer portion of ilie aperture. Thus, while all ffcequendes of the wput 
ooatiibute to the «hallo«wlip portion of the migrated image, only the penmssible low ft©- 
quendes of ae input data contribute to imaging the steepest dips. ^ ^ ^- 

Using a simple realiration of a ftequency-sdeetivB Kirdihoff migration operator, this 
technique is fllustfated on a synthetic data set iovolvtag greater than vertical dips. 

Introduction - 

For Integral migration methods, a major advantage of working in the time domain 
rather than tha frequency domain is speed, and a major disadvantage is inflentattV 
of the migration operator. In the frequency domain, the migration operator can be 
changed for each frequency, allowing for an optnnaJ combtnation 8*«P-diP 
imaging accuracy and «dection of unwanted data. For example, the combinattM^ of 
steep dips, high frequencies and large CMP spacing presents no serious processing 
problems for integral migration performed in the fcequenqr ^^i^^^^^^ 
frequencies, the migration aperture can be unlimited; for the highest frequenaea, the 

^ Received January 1991, revision accepted February 1992. Air iaim xix^k 

» Amooo Production Company. Reseen* Center. P.O. Box 3385. Tulsa. OK 74102. USX 
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aperture must be limited by the need to exclude aliased data from tbe migration at 
steep dips. By contrast, time-domain integral migration is performed only once on a 
given section (not once per frequency), but the aperture must be limited to the 
worst-case fiequency-domain nngration aperture, namely the one used at the 
highest frequencies. Usodting tbe aperture in this t^ay can result in tbe loss of 4 
significant portion of tbe low«freqnency steep-dip information in the migration. Not 
limiting tbe aperture can result in operator aliasing, which appears in the form 6t 
arti&cts caused by migrating high-frequency data past their maximum allowable 
dip. 

The most conmionly applied solution to this problem, short of migrating in the 
frequency domain, is to interpolate the input seismic data to a finer spatial grid' 
before migration. This results in an input section with smaller CMP spacing than 
the original section, and it increases the aliasing frequency for tihe migration accord- 
ing to the relation 

/«axlsinfl| 1 /J) 
v{0) ^4Ax' 

where is the m^*^"^"^ frequency retained in the data before migration, 9 is the 
angle of emergence of a reflection event, t<0) is the upper-sur£ace acoustic velocity 
and Ax is the CMP spadng. This inequality is a statement of tbe Nyquist criterion 
for sampling, on the discrete grid of CMP locations, a continuous wave-field areated 
by exploding reflectors in the earth: for each frequency/, unaliased transverse wave- 
numbers in the sampled wavefield (29^sin 8 divided by half-velocity at the earth's 
suiface) are limited in magnitude by 2n/2doc. To ensure that a KirchhofF migration 
handles only unaliased input data, the Nyquist criterion must hold simultaneously ' 
for all temporal frequencies, so tixat/^ appears in (1)- Attempting to migrate reflec- 
tion data witii temporal frequencies greater than/^ out to steep dips can result in 
artifacts which, in the worst case, might be interpreted as geological stnictura 
Decreasing Ax. increases the range of dips to which a giv» frequency component 
can be migrated; alternatively, it increases the range of frequencies which may con- 
tribute to a given dip in the image. Interpolating seismic data to decrease the CMP 
spacing, for the purpose of improving the migration operator response^ requires an 
accurate interpolation routine (to avoid creating additional aliased dati), and it 
increases the size of the data volume to be migrated^ hopefully without creating new 
information. The cost of the interpolation is usually much smaller than the migra- 
tion cost, and the added data volmne increases the migration cost, though not 
usually to the point of being more expensive than a frequency^-domain migration 
performed on the uninterpolated section. On a mainframe computer, the process of 
spatially interpolating a CMP section and then migrating using a time-domain inte- 
gral method is often (but not always) an acceptable way to obtain a high*quality 
migrated section in 2D with a minimum of migmtion axti&cts. 

However, in 3D, multiplicative factors of two in computation and data storage 
costs can critically affect the choice of processing parameters for a post-stack migra- 
tion. Target structures in 3D usually have their most rapid lateral geological varia- 
tion in a preferred direction; this is usually the direction of choice for a single 2D 
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gratipn at seismic line, and it usually defines the in-line, as opposed to cross-line, direction for 

once on a a 3D survey. Current practice for 3D acquisition and processing calls for the grid 

ed to the ' spadxkg in the cross-Iine direction to exceed the grid spacing in the in-line direction, 

ad at the often by a factor greater than two. Premigration processing may be unaffected by 

loss of a tius difference in spadngs, but 3D migration can be affected draniatically* Therefore^ 

^tion. Not ^ attempt to make the migration grid isotropic and to minimize migration arti- 

le form of . . interpolation is often performed in the cross-line direction after stack and 

allowable I before migration. 

The purpose of this note is to examine an alternative to data interpolation which 

ting in the can mitigate the problems of migrating seismic data with a large CMP spacing. This 

^atial grid method is applicable to time^domain integral migration, (In frequency-domain inte- 

udng than 8^ migration, it is performed automatically). Because it does not require data inter- 

^n accord- polation^ computation and data storage costs are reduced Thus, the proposed 
method has the potential to allow 3D migrations on large data sets without com- 
- promising migration quality. 

(1) 

MODIFYlliTG THE KlRCHHOPP MIGRATION OPERATOR 

ic velocity 

it criterion Frequency-domain migration avoids operator aliasing by restricting the migration 
sid created aperture when migrating the highest fr^uendes. In effect^ this procedure ignores the 
erse wave- portion of the input data with high transverse wavenumbers at high finequendes. A 
the earth's filter could be designed to perform this operation aplidtly in the frequency- 
'migration wavenumber domain before a time-domain migration, but the cost of the double 
iltaneously ' Fourier transform and its inverse can be comparable to the cost of the spatial inter- 
rate reflec- polation. Instead, the effect of the firequency-domatn migration operator selection 
n result in can be simulated very simply. Given an aperture suitable for migrating out to a 
structure. 1 specified maxnnum dip, the frequendes in the data whidi migrate out to tibe fiurthest 
.^mponent offsets of the migration aperture axe restricted and all frequencies whicb migrate out 
I may con- to the shorter o&ets axe left intact Thus, the inner portion of a diflEraction curve 
* the CMP (Fig. 1) refers to the smallest dips. None of the fi^uendes in the seismic data are 
equires an aliased for such small dips, so we can retain all the &equendes when migrating over 
taX and it that portion of the dijBpraction curve. On the other hand, the outer portion of the 
eating new ^ diflfraction curve rtfers to the steepest geological dips, for which the hi^best £re- 
the migra- quendes present in the r^ection events might be aliased When migrating over this 
hough not portion of the diffraction curve, oxxly the lowest, unaliased firequendes may be 
migration retained. This requires a second copy of the input traces, containing no alksed 
process of firequendes. (The aliashig firequency can be computed firom (1\ using the an^ of 
imain inte- -miergence of the ray which attains the steepest dip.) The third portion of the difi^ac- 
igh-quality tion curve, referring to intermediate dips, should be migrated using a third copy of 
the input traces, whose Ugh-firequenpy content is less than that of the original traces 
^ta storage but greater than that of the second copy. If the migration proceeds by smearing out 
ack migra* input traces onto the output section, rather than by summing along diffraction 
,ical varia- curves, the extra trace copies can be created, used and discarded along with each 
single 2D origmal input trace at almost no extra processing or storage cost (bi 3D migration. 
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OFFSET 

Fic., 1. Raypaihs through a depdi*varying medium, defining the portions of the diffiraclion 
carve for fhe various filtered versions of an input trace. The inner raypath, with the smaller 
take-off angle from the v^eal» bounds the segment of the diftaction curve for which the 
input trace needs no reduction in high«£:equency content. The outer raypath bounds the 
segment of the dxfiGraction curve for which the input trace needs the most reduction in hi^* 
frequency content 
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the basic unit of input data might be a block of several traces. This process then 
requires two extra blocks for a single block of original input traces, resulting in a 
non-triyial extra storage requirement) 

Id its crudest form, the modified migration operator requires the creation of two 
or more extra copies of each input trace. The second copy has all the aliased fre- 
qucndes removed, requiring a more severe high-cut filter than the one applied to the 
original traces. Ihe third copy can be an average of the first two, retaining the 
highest frequend^ at reduced amplitude. This approadbt yields a particular type of 
spatially-vaiying filter on the iiiput traces which is analogous to some applications 
of time-varying temporal filters. That is, a time-vaiymg filtered trace can be con- 
structed by a weighted average, at every time sample, of several dififerent band-pass- 
filiered versions of that trace; usually, lower-fi:equency versions are weighted more 
heavily at later times. Similarly, a spatially filtered set of input traces can be fonned 
as weighted averages of a few band-pass-filtered versions of a single trace, with the 
lower^firequency versions being weii^tcd more heavily at the largest offsets of the 
migration curvei Because the migration operator chooses a particular interpolated 
trace for mapping onto each oflSset, the procedure can be regarded as a migration 
operator interpolation. 
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The amount of aperture swept out by each of the (say, three) copies of an input 
trace is deteAnined by the take-off angle of the ray sent from the trace location on 
the upper surface to the hnaging location at dq)th. Figure 1 illustrates this. Ray- 
paths with take-off angles obeying (1) with/„^ of the original trace define the 
portion of the aperture swept out by the original trace, which has the greatest high- 
frequency content; this portion is lin&ited to the small oflbets. Sinuiar]^^, the outer 
portion of the diflBraction curve, determined by raypaths with the largest talce-off 
angles specified by the migration program, requires a reduction offg^ in order to 
obey (1). The second copy of the input trace, wHidh sweeps out the outer portion, 
will thus be the most severely high-cut version. The third copy of the input trace, 
which is an average of the first two, sweeps out the intermediate portion of the 
apwture. This copy of the trace can have a non-zero content of aliased frequencies, 
which can cause some problems when it sweeps out the &rthest offsets of the inter- 
mediate portion of the migration aperture. Therefore, it may be desirable either to 
reduce slightly the maxunum ofl&et sw^t out by this trace copy or to perfomi a 
more sophisticated interpolation whw sweeping out the intermediate offsets. 

Yihnaz (1987, p. 324) mentions an alternative operator modification which is 
sometimes used For the outer portion of the difi&action curve, more than one 
sample per trace is included for the summation which approximates the migration 
integral Hale (1991) provides specific rules for applying this modification to integial 
dq) moveout This approach will have the same efi^ as the method proposed here 
(that is, it will pezfoixn the same spatial filtering operation on the input -traces, 
resulting in the reduction of arti£acts due to migrating aliased data). However, it is 
clearly more expensive to implement, since it adds operations to the most time- 
consuming part of the migration. 



An Example 

Figure 2 shows a finite-difference zero-ofiset synthetic seismogram from a salt dome 
with an overhung face. Acoustic velocities in the outlying sediments are typical for 
the Gulf of Mexico, allowing reflection energy from the overhung portion of the salt 
flank to turn upwards and be recorded. CMP spadang is 33 m and the maximum 
frequency present in the data is approximately 30 Hz. Figure 3 shows a tumed^ray 
KirchhoJBf migration of the input data, where no effort has been made to suppitss 
the aliased frequencies at the steepest dips before or during migration. The migra- 
tion was performed in depth and stretched back in tune for display purpose. The 
steeply-dippmg salt flank, including the overhung portion with a dip of 110% has 
been well imaged, but the top part of the section is obscured by artifects resulting 
from migrating the aliased frequencies. In feet, high-frequency data from the flat 
reflectors have been migrated to steep dips which obscure the image of the sedi* 
ments. The steep-dip arti&cts appearing inside and immediatdy above the salt dome 
result from migrating the myriad of steep-dip multiple reflection and reverberation 
events evident on Fig. 2; the modified migration operator is not designed to attack 
thm. Figure 4 shows the modified Kirchhoff migration of the data, where the 
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Rg. 1 Finite-diflFerOTce2Brcw)fl&etsyndiedcse^ 
The steepest dip is 110'*. 
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FlG- 4. Ttimed-ray Kirchhoflf migration, using a ftequenpy-sdcctive operator. The artiiaocs 
multing from migrating aliased data have been $appzt$$ed. 

nugration operator has been modified using the crude procedure described in the 
px0Viou$ section. Figures 3 and 4 are scaled identically. As desired, the migration 
artifacts arising from migrating aliased frequencies have almost completely disap' 
peared. In adcUtion, fliere is no visible loss of high-irequency content on the fiat 
reflectors relative to Fig. 3. 



A technique has been proposed to eliminate one of the major disadvantages of 
time*domain Kircshhoff migration, namely the infleidbaity of the inigration operator 
and its resultant inability to image steep dips while induding high frequencies. This 
technique is applicable both to time migration, which implicitly uses straigjit rays in 
its diffraction curve calculations, and to depth migration, which uses movo sophisti- 
cated ray tracing. It is also applicable to integral dip moveout, where operator alia- 
ang is a serious drawback. A major application should be in 3D migration, dither 
single-pass or two-pass, where data interpolation presently serves as a fairly expen- 
sive and not entirely satisfactory alternative to migration operator aliasing. 



Hale» D. 1991. A nonaliased int^al method for dip moveout Geophysics 56, 795-805. 
YltMAZ, 0. 1987. Seismk Data Processing. Sodety of Exploration Oeopbysidsts, Tuba. 
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# Time Domain 16-degree Diffraction Movie 

# Star: w=p{t ,2) y=p(t ,z+l) ^ " 

# Star: u=p(t+l,z) v=p(t+l,z+l) 
realp(36,96),u(36),w{36),v(38Xy(36),e(36),f(36Xd(36),z(96),alfa.beta ^ 
integer ix^,iz,nz,it,nt,kbyte 

nx =s 36; nz = 96; nt = 96; kbyte=l 
alfa = .125 # v*dz*dt/(8*dx*dx) 

beta « .140 # accurate x derivative parameter; simplest case b=0. 

open(3,file= *plot40',status~ 'new *,accesB= 'direct ',fonn='miformatted',recl= 1) 
do iz=l,n2; do ix=l,nx; p(ix«iz} = 0. # dear space 

do iz=nz/5,nz,nz/4 5^ Set up initial model 

do it=l,15 # of 4 band limited 

do ix=l,4 # "point" scatterers. 

p(ix,it-J-iz) = (5.-ix)*(8-it)*exp(-.l*(it-8)**2) 
apb = alfa+beta; amb = aJf a-beta # tridiagonal coefficients 

diag = l.+2.*amb; offdi -amb 

do iz=nz.2,-2 { # ciimb up m steps of 2 z-levels 

doi=l,nz; z(i)=0.; 2(iz)==l. # Pointer to current a-level 
write(3,rec^kbyte)(z(ili=l,nz),((p(ix,i),i=^l.nz),ii=l,nx) 
kbyte = kbyte + nx*nz*4 + nz*4 
do ix— l,nx 

{ «(ix) = p(ix,iz-l); v(fac) = u(ix) } 

do it=iz,nt { 

do ix=l,nx #update the differencing star 

{ w(ix) = u(ix); y(ix) = v(ix); v(bc) = p(ix,it) ) 
dd = (l.-apb)*(v(l)+w(l))+apb*(v(2)+w(2)) 
d(l) = dd-diag*y(l)-offdi*(y(l)+y(2)) 
do ix=2,nx-l { 

dd = (l.-2.*apb)*(v{ix)+w(ix)) 
dd = dd + apb*(v(ix-l)+wfix-l)+v(ix+l)+w(ix+l)) 
d(ix) = dd-dia€*y(ix>-offdi*(y(ix-l)+y(ix+l)) } 
dd «= (l--apb)*(v(nx)+w(nx))+apb*(v(nx-l)+w(nx-l)) 
d(nx) = dd-dia«*y(nx)H3ffdi*(y(nx)+y(nx-l)) 
call rtris(nx,diag+offdi,offdi,diag,offdi,diag+offdi,d,u,e,f) 
do ix— 1^ 

J p(ixjt) = u(ix) 

do i=l,nz; z(i)=0.; z(l)=l. 

write(3^ec=kbyte) (z(i),i=l,n2),((p(ix,i),i=l,nz),ix=l.nx) 
stop; end 

Time-domain difiraction movie program. (Clayton, Gonzalez, 

JF C, Hale) 



Figure 2 shows the last frame in the movie produced by the test pro- 
gram. Exercise 1 suggests minor changes to the program of figure 1 to con- 
vert it from diffraction to migration. As modified, the program is essentially 
the original wave equation migration program introduced by Johnson and 
Claerbout [1971] and Doherty and Claerbout.[l972]. 
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